Crowding

Fitting

###Parameter Analyses###
rm(list=ls())
library(ggplot2)
library(mixRSVP)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(magrittr)
library(BayesFactor)

## Loading required package: coda

## Loading required package: Matrix

## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

participantPlots = TRUE #plot histograms with density?

inclusionBF <- function(priorProbs, variable){
  
  ###https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp###
  
  
  if(typeof(priorProbs) == 'S4') priorProbs <- as.vector(priorProbs)
  
  
  theseNames <- names(priorProbs)
  nProbs <- 1:length(priorProbs)
  variableMatches <- grep(variable, theseNames)
  
  if(grepl(':', variable)){
    subordinateVariables <- variable %>% strsplit(':') %>% unlist()

    thisRegex <- paste0(subordinateVariables,collapse = '.*\\+.*')
    
    subordinateEffects <- grep(thisRegex, theseNames, perl = T)
    subordinateEffects <- subordinateEffects[!subordinateEffects %in% variableMatches]
    
    
    sum(priorProbs[variableMatches])/sum(priorProbs[subordinateEffects])
  } else {
    interactionMatches <- grep(paste0(variable,'(?=:)|(?<=:)',variable), theseNames, perl = T)
    
    variableMainEffects <- variableMatches[!variableMatches %in% interactionMatches]
    
    
    otherMainEffects <- nProbs[!nProbs %in% c(variableMainEffects,interactionMatches)]
    
    
    sum(priorProbs[variableMainEffects])/sum(priorProbs[otherMainEffects])
  }
}



allErrors <- read.table('../allErrors.txt', sep='\t', stringsAsFactors = F, header = T)

nReps <- 100

bounds <- parameterBounds()

bounds['precision','upper'] <- 3

runAnyway <- FALSE #If TRUE, fit models regardless of the presence of a parameter file. 
plots <- FALSE

nParamFiles <- length(list.files(path = '../../modelOutput/',pattern ='parameterEstimates.*\\.csv',full.names = T)) #How many parameter DFs are saved?

if(nParamFiles>0){
  paramFiles <- list.files(path = '../../modelOutput',pattern ='parameterEstimates.*\\.csv',full.names = T) #What are the saved param files?
  
  splits <- strsplit(paramFiles, 'Estimates|(?<=[0-9][0-9])_|\\.csv',perl = T) #split out the date stamps from the param files names
  
  theseDates <- lapply(splits,FUN =function(x) paste(x[2],x[3])) %>% 
    unlist %>% 
    as.POSIXct(.,format='%d-%m-%Y %H-%M-%S') #Format them as POSIXct (datetime) 
  
  whichMaxDate <- which(theseDates == max(theseDates)) #which is the newest?
  params <- read.csv(paramFiles[whichMaxDate]) #load the newest
  
  #Empty space in the param DF for the un-modelled participants
  notModelled <- allErrors %>% filter(., !ID %in% params$ID) %>% pull(ID) %>% unique
  
  unModelledRows <- expand.grid(
    ID = factor(notModelled),
    crowded = factor(unique(allErrors$crowded)),
    ring = factor(unique(allErrors$ring)),
    efficacy = 999,
    latency = 999,
    precision = 999,
    val = 999,
    valGuessing = 999,
    pLRtest = 999,
    stringsAsFactors = F
  )
  
  params <- rbind(params, unModelledRows)
  
} else {
  notModelled <- allErrors %>% pull(ID) %>% unique
  
  params <- expand.grid(
    ID = factor(notModelled),
    crowded = factor(unique(allErrors$crowded)),
    ring = factor(unique(allErrors$ring)),
    efficacy = 999,
    latency = 999,
    precision = 999,
    val = 999,
    valGuessing = 999,
    pLRtest = 999,
    stringsAsFactors = F
  )
}  

if(length(notModelled)>0 & runAnyway){
  for(thisParticipant in notModelled){
    for(thisCondition in unique(allErrors$crowded)){
      for(thisRing in unique(allErrors$ring)){
        #print(paste0('Participant: ', thisParticipant, '. Ring: ', thisRing, '. Condition: ', thisCondition))
        
        try({
          invisible( #invisible combined with capture output hides those print 'error' statements hidden deep in mixRSVP::
            capture.output(
              theseParams <- allErrors %>% filter(., crowded == thisCondition, ring == thisRing, ID == thisParticipant) %>% analyzeOneCondition(., 24, bounds, nReps)
            )
          )
        }, silent = T) #Doing everything I can to silence those messages
      
        if(theseParams$pLRtest<.05){
          params %<>%
            mutate(efficacy=replace(efficacy, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$efficacy)) %>%
            as.data.frame()
          
          params %<>%
            mutate(latency=replace(latency, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$latency)) %>%
            as.data.frame()
          
          params %<>%
            mutate(precision=replace(precision, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$precision)) %>%
            as.data.frame()
        } else {
          params %<>%
            mutate(efficacy=replace(efficacy, ID == thisParticipant & crowded == thisCondition & ring == thisRing, 0)) %>%
            as.data.frame()
          
          params %<>%
            mutate(latency=replace(latency, ID == thisParticipant & crowded == thisCondition & ring == thisRing, NaN)) %>%
            as.data.frame()
          
          params %<>%
            mutate(precision=replace(precision, ID == thisParticipant & crowded == thisCondition & ring == thisRing, NaN)) %>%
            as.data.frame()
        }
        
        params %<>%
          mutate(val=replace(val, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$val)) %>%
          as.data.frame()
        
        params %<>%
          mutate(valGuessing=replace(valGuessing, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$valGuessing)) %>%
          as.data.frame()
        
        params %<>%
          mutate(pLRtest=replace(pLRtest, ID == thisParticipant & crowded == thisCondition & ring == thisRing, theseParams$pLRtest)) %>%
          as.data.frame()
      }
    }
  }
  write.csv(params, paste0('modelOutput/parameterEstimates',format(Sys.time(), "%d-%m-%Y_%H-%M-%S"),'.csv'),row.names = F)
}

paramsForAnalysis <- params %>% filter(efficacy>.1 & ID != 'CH')
paramsForAnalysis %<>% filter(!latency >= bounds[2,2] & !latency <= bounds[2,1])
paramsForAnalysis %<>% filter(!precision >= bounds[3,2] & !precision <= bounds[3,1])

paramsForAnalysis$ring %<>% as.factor
paramsForAnalysis %<>% mutate(IDChar=ID) %>% mutate(ID = as.numeric(ID)) %>% mutate(ID=as.factor(ID)) %>% as.data.frame() #For some reason the string 'AC5SONA' is interpreted as a nul by anovaBF, which uses base64

Efficacy Analyses

efficacyBF <- anovaBF(efficacy ~ ring * crowded + ID, 
                      data=paramsForAnalysis,
                      whichRandom = 'ID'
                      ) 

efficacyPriorProbs <- efficacyBF %>% newPriorOdds() %>% `*`(efficacyBF) %>% as.BFprobability() %>% as.vector() #extract P(M|Data) from BF object

efficacyInclusionBFs <- rep(-999, times = 3)

names(efficacyInclusionBFs) <- c('crowded','ring','crowded:ring')

for(name in names(efficacyInclusionBFs)){ #calculate inclusion BFs for main effects and interactions
  thisInclusionBF <- inclusionBF(efficacyPriorProbs, name)
  efficacyInclusionBFs[name] <- thisInclusionBF
}

knitr::kable(efficacyInclusionBFs)

	x
crowded	0.4620900
ring	135.8714676
crowded:ring	0.2057757

#Only evidence for an effect of ring
ggplot(paramsForAnalysis, aes(x=crowded, y = efficacy))+
  geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
  geom_point(aes(group = factor(ring)), position = position_dodge(.9))+
  stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = efficacy))+
  geom_point(aes(colour = ID), position = position_dodge())+
  geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
  #stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  #stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
  lims(y=c(0,1))+
  facet_wrap(~crowded)

## Warning: Width not defined. Set with `position_dodge(width = ?)`

## Warning: Width not defined. Set with `position_dodge(width = ?)`

Latency analysis

latencyBF <- anovaBF(latency ~ ring * crowded + ID, 
                      data=paramsForAnalysis,
                      whichRandom = 'ID'
)

latencyPriorProbs <- latencyBF %>% newPriorOdds() %>% `*`(latencyBF) %>% as.BFprobability() %>% as.vector() #that nested newPriorOdds() function call is the shittiest hack. Prior odds are 1 by default. I shouldn't need it, but S4 is a mystery to me

latencyInclusionBFs <- rep(-999, times = 3)

names(latencyInclusionBFs) <- c('crowded','ring','crowded:ring')

for(name in names(latencyInclusionBFs)){
  thisInclusionBF <- inclusionBF(latencyPriorProbs, name)
  latencyInclusionBFs[name] <- thisInclusionBF
}

knitr::kable(latencyInclusionBFs)

	x
crowded	2.2341273
ring	0.4563018
crowded:ring	0.5005254

ggplot(paramsForAnalysis, aes(x=crowded, y = latency))+
  geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
  geom_jitter(aes(group = factor(ring)), position = position_dodge(.9))+
  stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = latency))+
  geom_point(aes(colour = ID), position = position_dodge())+
  geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
  #stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  #stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
  facet_wrap(~crowded)

## Warning: Width not defined. Set with `position_dodge(width = ?)`

## Warning: Width not defined. Set with `position_dodge(width = ?)`

Precision Analysis

precisionBF <- anovaBF(precision ~ ring * crowded + ID, 
                     data=paramsForAnalysis,
                     whichRandom = 'ID'
)

precisionPriorProbs <- precisionBF %>% newPriorOdds() %>% `*`(precisionBF) %>% as.BFprobability() %>% as.vector() #that nested newPriorOdds() function call is the shittiest hack. Prior odds are 1 by default. I shouldn't need it, but S4 is a mystery to me

precisionInclusionBFs <- rep(-999, times = 3)

names(precisionInclusionBFs) <- c('crowded','ring','crowded:ring')

for(name in names(precisionInclusionBFs)){
  thisInclusionBF <- inclusionBF(precisionPriorProbs, name)
  precisionInclusionBFs[name] <- thisInclusionBF
}

knitr::kable(precisionInclusionBFs)

	x
crowded	0.3971054
ring	2192.4479671
crowded:ring	0.2211942

ggplot(paramsForAnalysis, aes(x=crowded, y = precision))+
  geom_violin(aes(fill = factor(ring)), position = position_dodge(.9))+
  geom_jitter(aes(group = factor(ring)), position = position_dodge(.9))+
  stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))

ggplot(paramsForAnalysis, aes(x=factor(ring), y = precision))+
  geom_point(aes(colour = ID), position = position_dodge())+
  geom_line(aes(group = ID, colour = ID), position = position_dodge(), alpha = .7)+
  #stat_summary(geom = 'point', aes(group = factor(ring)),fun.y = mean, position = position_dodge(.9))+
  #stat_summary(geom= 'errorbar', aes(group = factor(ring)), fun.data = mean_se, position = position_dodge(.9))+
  facet_wrap(~crowded)

## Warning: Width not defined. Set with `position_dodge(width = ?)`

## Warning: Width not defined. Set with `position_dodge(width = ?)`

Plots with Density

paramsForAnalysis %<>% mutate(ID = IDChar) %>% as.data.frame #Convert the problem name back for plotting

if(participantPlots){
  for(thisParticipant in unique(paramsForAnalysis$ID)){
    theseErrors <- allErrors %>% filter(ID==thisParticipant)
    theseDensities <- data.frame(density = numeric(0),SPE = numeric(0), crowded = character(0), ring = numeric(0), stringsAsFactors = F)
    theseConditions <- paramsForAnalysis %>% filter(ID == thisParticipant) %>% pull(crowded) %>% unique
    for(thisCondition in theseConditions){
      theseRings <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition) %>% pull(ring) %>% unique
      for(thisRing in theseRings){
        thisEfficacy <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(efficacy)
        thisLatency <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(latency)
        thisPrecision <- paramsForAnalysis %>% filter(ID == thisParticipant & crowded == thisCondition & ring == thisRing) %>% pull(precision)
        
        thisConditionErrors <- theseErrors %>% filter(crowded == thisCondition & ring == thisRing)
        # print(paste0('Participant: ', thisParticipant, '. Ring: ', thisRing, '. Condition: ', thisCondition,'. N = ', nrow(theseErrors)))
        minError <- thisConditionErrors %>% pull(SPE) %>% min
        maxError <- thisConditionErrors %>% pull(SPE) %>% max
        thisRange <- seq(minError,maxError,.1)
        
        thisConditionDensities <- data.frame(SPE = thisRange, density = dnorm(thisRange, thisLatency, thisPrecision), crowded = thisCondition, ring = thisRing, stringsAsFactors = F)
        theseDensities <- rbind(theseDensities, thisConditionDensities)
        # if(any(is.nan(theseDensities$density))){
        #   print(theseDensities)
        # }
        
        # thisFileName <- paste0('modelOutput/Plots/',thisCondition,'/',thisRing,'/',thisParticipant,'-',format(Sys.time(), "%d-%m-%Y_%H-%M-%S"),'.png')
        # ggsave(filename = thisFileName, thisPlot,width = 16, height = 9)
      }
    }
    thisPlot <- ggplot(theseErrors, aes(x=SPE))+
          geom_histogram(binwidth = 1)+
          geom_line(data = theseDensities, aes(x = SPE, y=density*60))+ #scale density to histogram with density * N * binwidth
          scale_y_continuous(sec.axis = sec_axis(~./60, name = 'Density'))+
          labs(y = 'Frequency', title = thisParticipant)+
          facet_wrap(~crowded*ring)
    show(thisPlot)
  }
}